#include "FEMSpace.h"
#include <typeinfo>
#include <math.h>
#include <ctime>
// #include "function1.h"
#include "function2.h"
using namespace std;
int main(int argc, char *argv[])
{
    // g++ -o main test_possion.cpp -std=c++11 -I /usr/include/eigen3/
    //  RectangleDomain* r = new RectangleDomain({{-1,-1},{1,-1},{1,1},{-1,1}});
    //  std::vector<Boundary<2> > B = r->boundary();
    vector<Point<2>> p = {{-1.0, -1.0}, {1.0, -1.0}, {1.0, 1.0}, {-1.0, 1.0}};
    Domain<2> *domain = new Rectangle2D(p);
    vector<int> seg({std::atoi(argv[1]), std::atoi(argv[2])});
    Mesh<2> *m = new Q2Mesh(domain, seg);
    TemplateElement<2> *e = new Quadrilateral_2_Element();
    Equation<2> *equ = new PossionEquation<2>();
    equ->SetRightHandsTermFunction(fe);
    equ->SetBoundaryConditionFunction({g1, g2, g3, g4});
    Possion_2D possionproblem(m, e, equ, domain);
    clock_t start, end;
    start = clock();
    possionproblem.setup_system();
    end = clock(); // 程序结束用时
    double endtime = (double)(end - start) / CLOCKS_PER_SEC;
    cout << "setup_system takes  " << endtime * 1000 << " ms" << endl; // ms为单位
    start = clock();
    possionproblem.assemble_system();
    end = clock(); // 程序结束用时
    endtime = (double)(end - start) / CLOCKS_PER_SEC;
    cout << "assemble_system takes  " << endtime * 1000 << " ms" << endl; // ms为单位
    start = clock();
    possionproblem.dealwithBoundaryCondition();
    end = clock(); // 程序结束用时
    endtime = (double)(end - start) / CLOCKS_PER_SEC;
    cout << "dealwithBoundaryCondition takes  " << endtime * 1000 << " ms" << endl; // ms为单位
    start = clock();
    possionproblem.solve();
    end = clock(); // 程序结束用时
    endtime = (double)(end - start) / CLOCKS_PER_SEC;
    cout << "solve takes  " << endtime * 1000 << " ms" << endl; // ms为单位

    Eigen::VectorXd Truesolution = Eigen::MatrixXd::Zero(m->getTotalNumDofs(), 1);
    Point<2> p0 = domain->getVertex(0);
    Point<2> p2 = domain->getVertex(2);
    double xLen = (p2[0] - p0[0]) / (1.0 * seg[0]) / 2;
    double yLen = (p2[1] - p0[1]) / (1.0 * seg[1]) / 2;
    for (int j = 0; j < 2 * seg[1] + 1; ++j)
    {
        for (int i = 0; i < 2 * seg[0] + 1; ++i)
        {
            // s[i+j*(seg[0]+1)] = (TrueSolve(p0[0]+i*xLen,p0[0]+j*yLen));//stupid!!!
            double point[2] = {p0[0] + i * xLen, p0[0] + j * yLen};
            Truesolution[i + j * (2 * seg[0] + 1)] = u(point); // stupid!!!
        }
    }
    Eigen::VectorXd residual = possionproblem.A() * possionproblem.getSolution() - possionproblem.Rhs();
    double residual_l2 = residual.squaredNorm();
    double residual_max = residual.lpNorm<Eigen::Infinity>();
    std::cout << "l2 residual is :" << residual_l2 << " ,residual max norm is :" << residual_max << std::endl;
    Eigen::VectorXd error = Truesolution - possionproblem.getSolution();
    double error_l2 = error.norm();
    double error_max = error.lpNorm<Eigen::Infinity>();
    std::cout << "Error l2 norm is :" << error_l2 << ", error max norm is :" << error_max << std::endl;
    if (std::atoi(argv[3]) == 0)
    {
        FILE *file1;
        file1 = fopen("output.txt", "w"); // “w”写文件覆盖原内容，不想覆盖则用“a”
        if (NULL == file1)
            printf("error");
        fprintf(file1, "%f\n%f\n%f\n%f\n", domain->getVertex(0)[0], domain->getVertex(1)[0], domain->getVertex(0)[1], domain->getVertex(3)[1]);
        fprintf(file1, "%d\n%d\n%d\n", m->getSegment()[0], m->getSegment()[1], 2);
        for (int i = 0; i < possionproblem.getSolution().size(); i++)
            fprintf(file1, "%.12f\n", possionproblem.getSolution()[i]); // 循环输出数据并写入
        fclose(file1);                                                  // 最后一定记得关闭文件
        std::cout << "Data has been output" << std::endl;
    }
    else if (std::atoi(argv[3]) == 1)
    {
        FILE *file2;
        file2 = fopen("err.csv", "a");
        if (NULL == file2)
            printf("error");
        fprintf(file2, "%d %d %.12f\n", std::atoi(argv[1]), std::atoi(argv[2]), error_l2);
        std::cout << "Data has been output" << std::endl;
    }
    else
        std::cout << "No output" << endl;
}
